from sympy import *
a=Symbol('a')
b=Symbol('b')
c=Symbol('c')
d=Symbol('d')
e=Symbol('e')
f=Symbol('f')
g=Symbol('g')
h=Symbol('h')
k=Symbol('k')

x=Symbol('x')
y=Symbol('y')



def basis(xh,a,b,c,d,e,f,g,h,k):
    return a+b*xh[0]+c*xh[1]+d*xh[0]*xh[1]+e*xh[0]**2+f*xh[1]**2+g*xh[0]**2*xh[1]+h*xh[0]*xh[1]**2+k*xh[0]**2*xh[1]**2


A=[
    [-1,-1],
    [0,-1],
    [1,-1],
    [1,0],
    [1,1],
    [0,1],
    [-1,1],
    [-1,0],
    [0,0]
]
print("p1")
f11=basis(A[0],a,b,c,d,e,f,g,h,k)-1
f12=basis(A[1],a,b,c,d,e,f,g,h,k)
f13=basis(A[2],a,b,c,d,e,f,g,h,k)
f14=basis(A[3],a,b,c,d,e,f,g,h,k)
f15=basis(A[4],a,b,c,d,e,f,g,h,k)
f16=basis(A[5],a,b,c,d,e,f,g,h,k)
f17=basis(A[6],a,b,c,d,e,f,g,h,k)
f18=basis(A[7],a,b,c,d,e,f,g,h,k)
f19=basis(A[8],a,b,c,d,e,f,g,h,k)

r=solve([f11,f12,f13,f14,f15,f16,f17,f18,f19],[a,b,c,d,e,f,g,h,k])
print(r)
f=r[a]+r[b]*x+r[c]*y+r[d]*x*y+r[e]*x**2+r[f]*y**2+r[g]*x**2*y+r[h]*x*y**2+r[k]*x**2*y**2
print(f)
print('df/dx')
print(f.diff(x,1,y,0))
print('df/dy')
print(f.diff(x,0,y,1))
print('d^2f/dxdy')
print(f.diff(x,1,y,1))
print('d^2f/dx^2')
print(f.diff(x,2,y,0))
print('d^2f/dy^2')
print(f.diff(x,0,y,2))
print('d^3f/dx^3')
print(f.diff(x,3,y,0))
print('d^3f/dy^3')
print(f.diff(x,0,y,3))
print('d^3f/dx^2dy')
print(f.diff(x,2,y,1))
print('d^3f/dxdy^2')
print(f.diff(x,1,y,2))
print('d^4f/dx^2dy^2')
print(f.diff(x,2,y,2))

print("\np2")

f21=basis(A[0],a,b,c,d,e,f,g,h,k)
f22=basis(A[1],a,b,c,d,e,f,g,h,k)-1
f23=basis(A[2],a,b,c,d,e,f,g,h,k)
f24=basis(A[3],a,b,c,d,e,f,g,h,k)
f25=basis(A[4],a,b,c,d,e,f,g,h,k)
f26=basis(A[5],a,b,c,d,e,f,g,h,k)
f27=basis(A[6],a,b,c,d,e,f,g,h,k)
f28=basis(A[7],a,b,c,d,e,f,g,h,k)
f29=basis(A[8],a,b,c,d,e,f,g,h,k)

r=solve([f21,f22,f23,f24,f25,f26,f27,f28,f29],[a,b,c,d,e,f,g,h,k])
print(r)
f=r[a]+r[b]*x+r[c]*y+r[d]*x*y+r[e]*x**2+r[f]*y**2+r[g]*x**2*y+r[h]*x*y**2+r[k]*x**2*y**2
print(f)
print('df/dx')
print(f.diff(x,1,y,0))
print('df/dy')
print(f.diff(x,0,y,1))
print('d^2f/dxdy')
print(f.diff(x,1,y,1))
print('d^2f/dx^2')
print(f.diff(x,2,y,0))
print('d^2f/dy^2')
print(f.diff(x,0,y,2))
print('d^3f/dx^3')
print(f.diff(x,3,y,0))
print('d^3f/dy^3')
print(f.diff(x,0,y,3))
print('d^3f/dx^2dy')
print(f.diff(x,2,y,1))
print('d^3f/dxdy^2')
print(f.diff(x,1,y,2))
print('d^4f/dx^2dy^2')
print(f.diff(x,2,y,2))

print("\np3")

f31=basis(A[0],a,b,c,d,e,f,g,h,k)
f32=basis(A[1],a,b,c,d,e,f,g,h,k)
f33=basis(A[2],a,b,c,d,e,f,g,h,k)-1
f34=basis(A[3],a,b,c,d,e,f,g,h,k)
f35=basis(A[4],a,b,c,d,e,f,g,h,k)
f36=basis(A[5],a,b,c,d,e,f,g,h,k)
f37=basis(A[6],a,b,c,d,e,f,g,h,k)
f38=basis(A[7],a,b,c,d,e,f,g,h,k)
f39=basis(A[8],a,b,c,d,e,f,g,h,k)

r=solve([f31,f32,f33,f34,f35,f36,f37,f38,f39],[a,b,c,d,e,f,g,h,k])
print(r)
f=r[a]+r[b]*x+r[c]*y+r[d]*x*y+r[e]*x**2+r[f]*y**2+r[g]*x**2*y+r[h]*x*y**2+r[k]*x**2*y**2
print(f)
print('df/dx')
print(f.diff(x,1,y,0))
print('df/dy')
print(f.diff(x,0,y,1))
print('d^2f/dxdy')
print(f.diff(x,1,y,1))
print('d^2f/dx^2')
print(f.diff(x,2,y,0))
print('d^2f/dy^2')
print(f.diff(x,0,y,2))
print('d^3f/dx^3')
print(f.diff(x,3,y,0))
print('d^3f/dy^3')
print(f.diff(x,0,y,3))
print('d^3f/dx^2dy')
print(f.diff(x,2,y,1))
print('d^3f/dxdy^2')
print(f.diff(x,1,y,2))
print('d^4f/dx^2dy^2')
print(f.diff(x,2,y,2))

print("\np4")

f41=basis(A[0],a,b,c,d,e,f,g,h,k)
f42=basis(A[1],a,b,c,d,e,f,g,h,k)
f43=basis(A[2],a,b,c,d,e,f,g,h,k)
f44=basis(A[3],a,b,c,d,e,f,g,h,k)-1
f45=basis(A[4],a,b,c,d,e,f,g,h,k)
f46=basis(A[5],a,b,c,d,e,f,g,h,k)
f47=basis(A[6],a,b,c,d,e,f,g,h,k)
f48=basis(A[7],a,b,c,d,e,f,g,h,k)
f49=basis(A[8],a,b,c,d,e,f,g,h,k)

r=solve([f41,f42,f43,f44,f45,f46,f47,f48,f49],[a,b,c,d,e,f,g,h,k])
print(r)
f=r[a]+r[b]*x+r[c]*y+r[d]*x*y+r[e]*x**2+r[f]*y**2+r[g]*x**2*y+r[h]*x*y**2+r[k]*x**2*y**2
print(f)
print('df/dx')
print(f.diff(x,1,y,0))
print('df/dy')
print(f.diff(x,0,y,1))
print('d^2f/dxdy')
print(f.diff(x,1,y,1))
print('d^2f/dx^2')
print(f.diff(x,2,y,0))
print('d^2f/dy^2')
print(f.diff(x,0,y,2))
print('d^3f/dx^3')
print(f.diff(x,3,y,0))
print('d^3f/dy^3')
print(f.diff(x,0,y,3))
print('d^3f/dx^2dy')
print(f.diff(x,2,y,1))
print('d^3f/dxdy^2')
print(f.diff(x,1,y,2))
print('d^4f/dx^2dy^2')
print(f.diff(x,2,y,2))

print("\np5")

f51=basis(A[0],a,b,c,d,e,f,g,h,k)
f52=basis(A[1],a,b,c,d,e,f,g,h,k)
f53=basis(A[2],a,b,c,d,e,f,g,h,k)
f54=basis(A[3],a,b,c,d,e,f,g,h,k)
f55=basis(A[4],a,b,c,d,e,f,g,h,k)-1
f56=basis(A[5],a,b,c,d,e,f,g,h,k)
f57=basis(A[6],a,b,c,d,e,f,g,h,k)
f58=basis(A[7],a,b,c,d,e,f,g,h,k)
f59=basis(A[8],a,b,c,d,e,f,g,h,k)

r=solve([f51,f52,f53,f54,f55,f56,f57,f58,f59],[a,b,c,d,e,f,g,h,k])
print(r)
f=r[a]+r[b]*x+r[c]*y+r[d]*x*y+r[e]*x**2+r[f]*y**2+r[g]*x**2*y+r[h]*x*y**2+r[k]*x**2*y**2
print(f)
print('df/dx')
print(f.diff(x,1,y,0))
print('df/dy')
print(f.diff(x,0,y,1))
print('d^2f/dxdy')
print(f.diff(x,1,y,1))
print('d^2f/dx^2')
print(f.diff(x,2,y,0))
print('d^2f/dy^2')
print(f.diff(x,0,y,2))
print('d^3f/dx^3')
print(f.diff(x,3,y,0))
print('d^3f/dy^3')
print(f.diff(x,0,y,3))
print('d^3f/dx^2dy')
print(f.diff(x,2,y,1))
print('d^3f/dxdy^2')
print(f.diff(x,1,y,2))
print('d^4f/dx^2dy^2')
print(f.diff(x,2,y,2))


print("\np6")

f61=basis(A[0],a,b,c,d,e,f,g,h,k)
f62=basis(A[1],a,b,c,d,e,f,g,h,k)
f63=basis(A[2],a,b,c,d,e,f,g,h,k)
f64=basis(A[3],a,b,c,d,e,f,g,h,k)
f65=basis(A[4],a,b,c,d,e,f,g,h,k)
f66=basis(A[5],a,b,c,d,e,f,g,h,k)-1
f67=basis(A[6],a,b,c,d,e,f,g,h,k)
f68=basis(A[7],a,b,c,d,e,f,g,h,k)
f69=basis(A[8],a,b,c,d,e,f,g,h,k)

r=solve([f61,f62,f63,f64,f65,f66,f67,f68,f69],[a,b,c,d,e,f,g,h,k])
print(r)
f=r[a]+r[b]*x+r[c]*y+r[d]*x*y+r[e]*x**2+r[f]*y**2+r[g]*x**2*y+r[h]*x*y**2+r[k]*x**2*y**2
print(f)
print('df/dx')
print(f.diff(x,1,y,0))
print('df/dy')
print(f.diff(x,0,y,1))
print('d^2f/dxdy')
print(f.diff(x,1,y,1))
print('d^2f/dx^2')
print(f.diff(x,2,y,0))
print('d^2f/dy^2')
print(f.diff(x,0,y,2))
print('d^3f/dx^3')
print(f.diff(x,3,y,0))
print('d^3f/dy^3')
print(f.diff(x,0,y,3))
print('d^3f/dx^2dy')
print(f.diff(x,2,y,1))
print('d^3f/dxdy^2')
print(f.diff(x,1,y,2))
print('d^4f/dx^2dy^2')
print(f.diff(x,2,y,2))


print("\np7")

f71=basis(A[0],a,b,c,d,e,f,g,h,k)
f72=basis(A[1],a,b,c,d,e,f,g,h,k)
f73=basis(A[2],a,b,c,d,e,f,g,h,k)
f74=basis(A[3],a,b,c,d,e,f,g,h,k)
f75=basis(A[4],a,b,c,d,e,f,g,h,k)
f76=basis(A[5],a,b,c,d,e,f,g,h,k)
f77=basis(A[6],a,b,c,d,e,f,g,h,k)-1
f78=basis(A[7],a,b,c,d,e,f,g,h,k)
f79=basis(A[8],a,b,c,d,e,f,g,h,k)

r=solve([f71,f72,f73,f74,f75,f76,f77,f78,f79],[a,b,c,d,e,f,g,h,k])
print(r)
f=r[a]+r[b]*x+r[c]*y+r[d]*x*y+r[e]*x**2+r[f]*y**2+r[g]*x**2*y+r[h]*x*y**2+r[k]*x**2*y**2
print(f)
print('df/dx')
print(f.diff(x,1,y,0))
print('df/dy')
print(f.diff(x,0,y,1))
print('d^2f/dxdy')
print(f.diff(x,1,y,1))
print('d^2f/dx^2')
print(f.diff(x,2,y,0))
print('d^2f/dy^2')
print(f.diff(x,0,y,2))
print('d^3f/dx^3')
print(f.diff(x,3,y,0))
print('d^3f/dy^3')
print(f.diff(x,0,y,3))
print('d^3f/dx^2dy')
print(f.diff(x,2,y,1))
print('d^3f/dxdy^2')
print(f.diff(x,1,y,2))
print('d^4f/dx^2dy^2')
print(f.diff(x,2,y,2))


print("\np8")

f81=basis(A[0],a,b,c,d,e,f,g,h,k)
f82=basis(A[1],a,b,c,d,e,f,g,h,k)
f83=basis(A[2],a,b,c,d,e,f,g,h,k)
f84=basis(A[3],a,b,c,d,e,f,g,h,k)
f85=basis(A[4],a,b,c,d,e,f,g,h,k)
f86=basis(A[5],a,b,c,d,e,f,g,h,k)
f87=basis(A[6],a,b,c,d,e,f,g,h,k)
f88=basis(A[7],a,b,c,d,e,f,g,h,k)-1
f89=basis(A[8],a,b,c,d,e,f,g,h,k)

r=solve([f81,f82,f83,f84,f85,f86,f87,f88,f89],[a,b,c,d,e,f,g,h,k])
print(r)
f=r[a]+r[b]*x+r[c]*y+r[d]*x*y+r[e]*x**2+r[f]*y**2+r[g]*x**2*y+r[h]*x*y**2+r[k]*x**2*y**2
print(f)
print('df/dx')
print(f.diff(x,1,y,0))
print('df/dy')
print(f.diff(x,0,y,1))
print('d^2f/dxdy')
print(f.diff(x,1,y,1))
print('d^2f/dx^2')
print(f.diff(x,2,y,0))
print('d^2f/dy^2')
print(f.diff(x,0,y,2))
print('d^3f/dx^3')
print(f.diff(x,3,y,0))
print('d^3f/dy^3')
print(f.diff(x,0,y,3))
print('d^3f/dx^2dy')
print(f.diff(x,2,y,1))
print('d^3f/dxdy^2')
print(f.diff(x,1,y,2))
print('d^4f/dx^2dy^2')
print(f.diff(x,2,y,2))



print("\np9")

f91=basis(A[0],a,b,c,d,e,f,g,h,k)
f92=basis(A[1],a,b,c,d,e,f,g,h,k)
f93=basis(A[2],a,b,c,d,e,f,g,h,k)
f94=basis(A[3],a,b,c,d,e,f,g,h,k)
f95=basis(A[4],a,b,c,d,e,f,g,h,k)
f96=basis(A[5],a,b,c,d,e,f,g,h,k)
f97=basis(A[6],a,b,c,d,e,f,g,h,k)
f98=basis(A[7],a,b,c,d,e,f,g,h,k)
f99=basis(A[8],a,b,c,d,e,f,g,h,k)-1

r=solve([f91,f92,f93,f94,f95,f96,f97,f98,f99],[a,b,c,d,e,f,g,h,k])
print(r)
f=r[a]+r[b]*x+r[c]*y+r[d]*x*y+r[e]*x**2+r[f]*y**2+r[g]*x**2*y+r[h]*x*y**2+r[k]*x**2*y**2
print(f)
print('df/dx')
print(f.diff(x,1,y,0))
print('df/dy')
print(f.diff(x,0,y,1))
print('d^2f/dxdy')
print(f.diff(x,1,y,1))
print('d^2f/dx^2')
print(f.diff(x,2,y,0))
print('d^2f/dy^2')
print(f.diff(x,0,y,2))
print('d^3f/dx^3')
print(f.diff(x,3,y,0))
print('d^3f/dy^3')
print(f.diff(x,0,y,3))
print('d^3f/dx^2dy')
print(f.diff(x,2,y,1))
print('d^3f/dxdy^2')
print(f.diff(x,1,y,2))
print('d^4f/dx^2dy^2')
print(f.diff(x,2,y,2))

print("\nanalytical solution for example 1")
ana=x*y*(1-x/2)*(1-y)*exp(x+y)
print('\ndy/dx')
print(ana.diff(x,1,y,0))
print('\ndy/dy')
print(ana.diff(x,0,y,1))
print('\nd^2y/dxdy')
print(ana.diff(x,1,y,1))

print("\nanalytical solution for example 2")
ana=exp(x+y)
print('\ndy/dx')
print(ana.diff(x,1,y,0))
print('\ndy/dy')
print(ana.diff(x,0,y,1))
print('\nd^2y/dxdy')
print(ana.diff(x,1,y,1))










